Set up workspace, set options, and load required packages.
rm(list=ls(all=TRUE))
knitr::opts_chunk$set(root.dir = "~/R",warning=FALSE, message=FALSE)
library(gridExtra)
library(multcomp)
library(emmeans)
library(ggplot2)
library(dplyr)
library(effects)
library(plyr)
library(lattice)
library(glmmTMB)
library(effects)
library(lsmeans)
library(MuMIn)
library(MASS)
library(car)
library(FSA)
library(lmerTest)
library(lme4)
library(blmeco)
library(ggsci)
library(coxme)
library(survival)
library(cowplot)
library(RColorBrewer)
library(tidyverse)
library(partitionBEFsp)
library(factoextra)
library(ggfortify)
library(cowplot)
library(survminer)
Survivorship is measured as binary on each individual spat in the experiment.
#load data
rearing<-read.csv("Data/GrowOut_Responses_PostGenotyping.csv", header=TRUE, sep=",", na.strings="NA")
rearing$Quadrant<-as.factor(rearing$Quadrant)
rearing$Diversity.Type<-relevel(rearing$Diversity.Type, ref="Individual")
#rearing$Total.Juveniles<-as.factor(rearing$Total.Juveniles)
#rearing$Fusion.Partners<-as.factor(rearing$Fusion.Partners)
#rearing$Richness<-as.factor(rearing$Richness)
Subset the final timepoint for additional analyses.
final<-rearing[which(rearing$Timepoint == "Final"),]
multiple<-rearing[which(rearing$Community == "Multiple"),]
multiple_final<-multiple[which(multiple$Timepoint == "Final"),]
Analyze with a logistic binomial regression model.
Analyze at the final timepoint.
#with final dataset
survmod2<-glmer(Survivorship ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family=binomial)
summary(survmod2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Survivorship ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners +
## (1 | Parent.Site/Colony) + (1 | Tank) + (1 | Slide)
## Data: final
##
## AIC BIC logLik deviance df.resid
## 224.9 255.5 -103.5 206.9 212
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.6175 -0.6293 0.3049 0.5167 1.4715
##
## Random effects:
## Groups Name Variance Std.Dev.
## Slide (Intercept) 3.186e-01 0.5644532
## Colony:Parent.Site (Intercept) 1.488e-02 0.1219734
## Parent.Site (Intercept) 1.634e-07 0.0004042
## Tank (Intercept) 1.709e-01 0.4133998
## Number of obs: 221, groups:
## Slide, 115; Colony:Parent.Site, 21; Parent.Site, 6; Tank, 5
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.450767 0.519055 2.795 0.00519 **
## Richness -0.007283 0.267024 -0.027 0.97824
## CommunityMultiple -1.847623 0.618356 -2.988 0.00281 **
## Fusion.Partners 3.036311 0.946317 3.209 0.00133 **
## Richness:Fusion.Partners -0.432505 0.319527 -1.354 0.17587
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Rchnss CmmntM Fsn.Pr
## Richness -0.503
## CmmntyMltpl -0.231 -0.556
## Fusn.Prtnrs -0.105 0.414 -0.505
## Rchnss:Fs.P 0.145 -0.474 0.405 -0.902
## convergence code: 0
## Model failed to converge with max|grad| = 0.00237872 (tol = 0.001, component 1)
Anova(survmod2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Survivorship
## Chisq Df Pr(>Chisq)
## Richness 0.5759 1 0.447925
## Community 8.9279 1 0.002808 **
## Fusion.Partners 21.1261 1 4.3e-06 ***
## Richness:Fusion.Partners 1.8322 1 0.175871
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#final dataset
summary(effect(c("Richness", "Fusion.Partners"), xlevels=4, survmod2))
##
## Richness*Fusion.Partners effect
## Fusion.Partners
## Richness 0 1 2 3
## 1 0.5201784 0.9361096 0.9949754 0.9996265
## 2 0.5183603 0.9041990 0.9880627 0.9986242
## 3 0.5165418 0.8587551 0.9719082 0.9949464
## 4 0.5147229 0.7966036 0.9353234 0.9816172
##
## Lower 95 Percent Confidence Limits
## Fusion.Partners
## Richness 0 1 2 3
## 1 0.3287517 0.8160149 0.9478899 0.9855863
## 2 0.3643003 0.7954607 0.9392540 0.9827277
## 3 0.3149904 0.7068390 0.8811552 0.9520535
## 4 0.2323719 0.5215122 0.6410654 0.7004323
##
## Upper 95 Percent Confidence Limits
## Fusion.Partners
## Richness 0 1 2 3
## 1 0.7058580 0.9797581 0.9995363 0.9999905
## 2 0.6690068 0.9581692 0.9977482 0.9998920
## 3 0.7128515 0.9387678 0.9938440 0.9994880
## 4 0.7879786 0.9336587 0.9915323 0.9991807
summary(effect("Fusion.Partners", xlevels=4, survmod2))
##
## Fusion.Partners effect
## Fusion.Partners
## 0 1 2 3
## 0.5186977 0.9110372 0.9898281 0.9989197
##
## Lower 95 Percent Confidence Limits
## Fusion.Partners
## 0 1 2 3
## 0.3642350 0.8027598 0.9426089 0.9840155
##
## Upper 95 Percent Confidence Limits
## Fusion.Partners
## 0 1 2 3
## 0.6696673 0.9626407 0.9982685 0.9999280
summary(effect("Community", xlevels=4, survmod2))
##
## Community effect
## Community
## Individual Multiple
## 0.9480743 0.7421160
##
## Lower 95 Percent Confidence Limits
## Community
## Individual Multiple
## 0.8373723 0.5918728
##
## Upper 95 Percent Confidence Limits
## Community
## Individual Multiple
## 0.9847894 0.8509763
dispersion_glmer(survmod2) #no overdispersion
## [1] 0.9357489
plot(residuals(survmod2)) + abline(h=0, lty=2) #Dispersed randomly
## integer(0)
Create plot for effect of community type.
eff.surv1 <- predictorEffect("Community", survmod2)
surv.effects1<-plot(eff.surv1,
lines=list(multiline=TRUE,
col=c("black"),
lty=1),
confint=list(style="bars", width=4, col="black"),
lwd=4,
axes=list(y=list(lim=c(0, 1))),
type="response",
ylab=expression(bold("Prob(Survivorship)")),
legend.position="right",
xlab=expression(bold("Number of Juveniles")),
main=""); surv.effects1
#lattice=list(key.args=list(space="right",
#border=FALSE,
#title=expression(bold("Genotypic \nRichness")),
#cex=1,
#cex.title=1)));surv.effects1
Create plots for fusion partner effects.
survmod2a<-glmer(Survivorship ~ Richness + Fusion.Partners + Richness:Fusion.Partners + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family=binomial, subset=c(Community=="Multiple"))
eff.surv2 <- predictorEffect("Fusion.Partners", survmod2a, xlevels=4, rug=TRUE)
surv.effects2<-plot(eff.surv2,
lines=list(multiline=TRUE,
col=c("lightblue","mediumblue", "blue4", "darkblue"),
lty=1),
confint=list(style="bands", alpha=0.1),
lwd=4,
axes=list(y=list(lim=c(0, 1))),
type="response",
ylab=expression(bold("Prob(Survivorship)")),
legend.position="right",
xlab=expression(bold("Fusion Partners")),
main="",
lattice=list(key.args=list(space="right",
border=FALSE,
title=expression(bold("Genotypic \nRichness")),
cex=1,
cex.title=1)));surv.effects2
survmod2b<-glmer(Survivorship ~ Fusion.Partners + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family=binomial, subset=c(Community=="Multiple"))
eff.surv2b <- predictorEffect("Fusion.Partners", survmod2b, xlevels=4, rug=TRUE)
surv.effects2b<-plot(eff.surv2b,
lines=list(multiline=TRUE,
col=c("black"),
lty=1),
confint=list(style="bands", alpha=0.1),
lwd=4,
axes=list(y=list(lim=c(0, 1))),
type="response",
ylab=expression(bold("Prob(Survivorship)")),
legend.position="right",
xlab=expression(bold("Fusion Partners")),
main="",
lattice=list(key.args=list(space="none",
border=FALSE,
title=expression(bold("")),
cex=1,
cex.title=1)));surv.effects2b
Generate mean values of survivorship for 1) number of juveniles and 2) within “multiple” juvenile groups, mean values for number of fusion partners.
meansurv <- ddply(final, c("Community"), summarise,
N = length(Survivorship[!is.na(Survivorship)]),
mean = mean(Survivorship, na.rm=TRUE),
sd = sd(Survivorship, na.rm=TRUE),
se = sd / sqrt(N),
max = max(Survivorship, na.rm=TRUE),
lower = mean-sd,
upper = mean+sd
)
meansurv
## Community N mean sd se max lower upper
## 1 Individual 58 0.7758621 0.4206554 0.05523477 1 0.3552066 1.196518
## 2 Multiple 163 0.6687117 0.4721270 0.03697984 1 0.1965847 1.140839
meansurv2 <- ddply(multiple_final, c("Fusion.Partners"), summarise,
N = length(Survivorship[!is.na(Survivorship)]),
mean = mean(Survivorship, na.rm=TRUE),
sd = sd(Survivorship, na.rm=TRUE),
se = sd / sqrt(N),
max = max(Survivorship, na.rm=TRUE),
lower = mean-sd,
upper = mean+sd
)
meansurv2
## Fusion.Partners N mean sd se max lower upper
## 1 0 75 0.3466667 0.4791133 0.05532324 1 -0.1324466 0.825780
## 2 1 48 0.9583333 0.2019409 0.02914766 1 0.7563924 1.160274
## 3 2 24 0.9166667 0.2823299 0.05763034 1 0.6343368 1.198997
## 4 3 16 0.9375000 0.2500000 0.06250000 1 0.6875000 1.187500
Analyze with a logistic binomial regression model.
Final timepoint dataset.
Analyze fusion as a product of genotypic richness within “multiple” juvenile groups.
#with final dataset
fusemod1<-glmer(Fusion ~ Richness + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family=binomial, subset=c(Community=="Multiple"))
summary(fusemod1)
Anova(fusemod1, type="II") # no effects on rate of fusion by the end point
dispersion_glmer(fusemod1) #no overdispersion
plot(residuals(fusemod1)) + abline(h=0, lty=2)
There is an effect of genotypic diversity on fusion rates. Plot the relationship.
eff.fuse1 <- predictorEffect("Richness", fusemod1, xlevels=4, rug=TRUE)
fuse.effects1<-plot(eff.fuse1,
lines=list(multiline=TRUE,
col=c("black"),
lty=1),
confint=list(style="bands", alpha=0.1),
lwd=4,
#axes=list(y=list(lim=c(0, 1))),
type="response",
ylab=expression(bold("Prob(Successful Fusion)")),
xlab=expression(bold("Richness")),
main="",
lattice=list(key.args=list(space="right",
border=FALSE,
title=expression(bold("Genotypic \nRichness")),
cex=1,
cex.title=1)));fuse.effects1
Generate a summary of proportion fused at the end of the grow-out period.
meanfusion <- ddply(multiple_final, c("Richness"), summarise,
N = length(Fusion[!is.na(Fusion)]),
mean = mean(Fusion, na.rm=TRUE),
sd = sd(Fusion, na.rm=TRUE),
se = sd / sqrt(N),
max = max(Fusion, na.rm=TRUE),
lower = mean-sd,
upper = mean+sd
)
meanfusion
## Richness N mean sd se max lower upper
## 1 1 59 0.5423729 0.5024778 0.06541704 1 0.03989509 1.0448507
## 2 2 52 0.7115385 0.4574670 0.06343925 1 0.25407150 1.1690054
## 3 3 28 0.1428571 0.3563483 0.06734350 1 -0.21349118 0.4992055
## 4 4 24 0.6250000 0.4945354 0.10094661 1 0.13046464 1.1195354
Growth is measured as polyp expansion in each spat by number of spat in each sibling type and number of fused individuals. Growth is measured in final dataset as all spat started with 1 polyp each.
Analyze polyp expansion with the final dataset and present means.
polypmod1<-glmer(Polyps ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family = poisson)
summary(polypmod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Polyps ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners +
## (1 | Parent.Site/Colony) + (1 | Tank) + (1 | Slide)
## Data: final
##
## AIC BIC logLik deviance df.resid
## 590.6 618.1 -286.3 572.6 148
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1240 -0.5469 -0.1021 0.4563 1.4273
##
## Random effects:
## Groups Name Variance Std.Dev.
## Slide (Intercept) 1.858e-09 4.311e-05
## Colony:Parent.Site (Intercept) 7.237e-02 2.690e-01
## Parent.Site (Intercept) 1.445e-03 3.802e-02
## Tank (Intercept) 2.789e-03 5.282e-02
## Number of obs: 157, groups:
## Slide, 92; Colony:Parent.Site, 21; Parent.Site, 6; Tank, 5
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.279359 0.147887 8.651 <2e-16 ***
## Richness -0.021274 0.091952 -0.231 0.817
## CommunityMultiple -0.095335 0.160455 -0.594 0.552
## Fusion.Partners 0.068683 0.143007 0.480 0.631
## Richness:Fusion.Partners -0.009974 0.052539 -0.190 0.849
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Rchnss CmmntM Fsn.Pr
## Richness -0.623
## CmmntyMltpl 0.044 -0.611
## Fusn.Prtnrs -0.291 0.579 -0.673
## Rchnss:Fs.P 0.384 -0.723 0.582 -0.912
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Anova(polypmod1, type="II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Polyps
## Chisq Df Pr(>Chisq)
## Richness 0.2843 1 0.5939
## Community 0.3530 1 0.5524
## Fusion.Partners 0.5587 1 0.4548
## Richness:Fusion.Partners 0.0360 1 0.8494
qqPlot(residuals(polypmod1))
## 389 393
## 73 74
meangrowth <- ddply(final, c("Community"), summarise,
N = length(Polyps[!is.na(Polyps)]),
mean = mean(Polyps, na.rm=TRUE),
sd = sd(Polyps, na.rm=TRUE),
se = sd / sqrt(N),
max = max(Polyps, na.rm=TRUE),
lower = mean-sd,
upper = mean+sd
)
meangrowth
## Community N mean sd se max lower upper
## 1 Individual 45 3.533333 1.778661 0.2651472 8 1.754672 5.311995
## 2 Multiple 112 3.357143 1.558800 0.1472928 7 1.798343 4.915943
Plot the model outputs.
eff.polyps1 <- predictorEffect("Fusion.Partners", polypmod1, xlevels=4, rug=TRUE)
polyp.effects1<-plot(eff.polyps1,
lines=list(multiline=TRUE,
col=c("lightblue","mediumblue", "blue4", "darkblue"),
lty=1),
confint=list(style="bands", alpha=0.1),
lwd=4,
#axes=list(y=list(lim=c(0, 1))),
type="response",
ylab=expression(bold("Total Polyp Expansion")),
legend.position="right",
xlab=expression(bold("Fusion Partners")),
main="",
lattice=list(key.args=list(space="right",
border=FALSE,
title=expression(bold("Genotypic \nRichness")),
cex=1,
cex.title=1)));polyp.effects1
Combine figures for Grow-Out Period.
grow_out_figs<-plot_grid(surv.effects1, surv.effects2b, labels = c("A", "B"), ncol=2, nrow=1, rel_heights= c(1,1), rel_widths = c(0.8,1), label_size = 20, label_y=1, align="h");grow_out_figs
ggsave(filename="Figures/grow_out_figs.png", plot=grow_out_figs, dpi=500, width=10, height=6, units="in")
Survivorship is measured as binary on each individual spat in the experiment.
#load data
stress<-read.csv("Data/Stress_Responses_PostGenotyping.csv", header=TRUE, sep=",", na.strings="NA")
stress$Diversity.Type<-relevel(stress$Diversity.Type, ref="Individual")
stress$Polyps<-as.numeric(stress$Polyps)
Subset the final timepoint, “multiple” juvenile datasets for additional analyses.
final_stress<-stress[which(stress$Timepoint == "End"),]
multiple_stress<-stress[which(stress$Community == "Multiple"),]
multiple_final_stress<-multiple_stress[which(multiple_stress$Timepoint == "End"),]
d26_stress<-stress[which(stress$Timepoint == "S2D11"),]
Analyze with a logistic binomial regression model.
Check for effect of community type (“individual” or “multiple” juveniles) on survivorship and display means.
#with full time dataset
survmod3b<-glmer(Survivorship ~ Day * Treatment * Community + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=stress, family=binomial, link=logit)
summary(survmod3b)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Survivorship ~ Day * Treatment * Community + (1 | Parent.Site/Colony) +
## (1 | Treatment:Tank) + (1 | Slide)
## Data: stress
##
## AIC BIC logLik deviance df.resid
## 449.0 508.9 -212.5 425.0 1076
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -110.079 0.001 0.028 0.118 5.533
##
## Random effects:
## Groups Name Variance Std.Dev.
## Slide (Intercept) 2.605e+00 1.614e+00
## Colony:Parent.Site (Intercept) 5.717e-01 7.561e-01
## Treatment:Tank (Intercept) 5.061e-09 7.114e-05
## Parent.Site (Intercept) 1.645e-08 1.282e-04
## Number of obs: 1088, groups:
## Slide, 90; Colony:Parent.Site, 20; Treatment:Tank, 9; Parent.Site, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.7660321 1.7130404 3.950 7.82e-05 ***
## Day -0.0835370 0.0592967 -1.409 0.158896
## TreatmentHigh 0.8459104 1.9259267 0.439 0.660500
## CommunityMultiple 3.4775795 2.4771976 1.404 0.160368
## Day:TreatmentHigh -0.2886216 0.0823811 -3.503 0.000459 ***
## Day:CommunityMultiple -0.1855821 0.0878809 -2.112 0.034708 *
## TreatmentHigh:CommunityMultiple 2.8355215 3.2070137 0.884 0.376608
## Day:TreatmentHigh:CommunityMultiple 0.0001848 0.1226640 0.002 0.998798
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Day TrtmnH CmmntM Dy:TrH Dy:CmM TrH:CM
## Day -0.835
## TreatmntHgh -0.700 0.666
## CmmntyMltpl -0.586 0.538 0.484
## Dy:TrtmntHg 0.397 -0.633 -0.837 -0.348
## Dy:CmmntyMl 0.507 -0.652 -0.450 -0.917 0.471
## TrtmntHg:CM 0.412 -0.395 -0.600 -0.760 0.504 0.702
## Dy:TrtmH:CM -0.295 0.436 0.563 0.641 -0.648 -0.713 -0.924
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Anova(survmod3b, type="II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Survivorship
## Chisq Df Pr(>Chisq)
## Day 65.5353 1 5.708e-16 ***
## Treatment 37.2722 1 1.027e-09 ***
## Community 3.3205 1 0.068421 .
## Day:Treatment 21.1660 1 4.212e-06 ***
## Day:Community 9.0584 1 0.002615 **
## Treatment:Community 5.3796 1 0.020374 *
## Day:Treatment:Community 0.0000 1 0.998798
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dispersion_glmer(survmod3b) #no overdispersion
## [1] 0.5599041
plot(residuals(survmod3b)) + abline(h=0, lty=2) #Dispersed randomly
## integer(0)
meanstresssurv1 <- ddply(final_stress, c("Community", "Treatment"), summarise,
N = length(Survivorship[!is.na(Survivorship)]),
mean = mean(Survivorship, na.rm=TRUE),
sd = sd(Survivorship, na.rm=TRUE),
se = sd / sqrt(N),
max = max(Survivorship, na.rm=TRUE),
lower = mean-sd,
upper = mean+sd
)
meanstresssurv1
## Community Treatment N mean sd se max lower
## 1 Individual Ambient 21 0.95238095 0.2182179 0.04761905 1 0.7341631
## 2 Individual High 24 0.00000000 0.0000000 0.00000000 0 0.0000000
## 3 Multiple Ambient 46 0.76086957 0.4312660 0.06358670 1 0.3296036
## 4 Multiple High 65 0.04615385 0.2114510 0.02622727 1 -0.1652972
## upper
## 1 1.1705988
## 2 0.0000000
## 3 1.1921355
## 4 0.2576049
Analyze within “multiple” groups as survival was significantly different between single and multiple juvenile slides.
#with full time dataset
survmod3<-glmer(Survivorship ~ Day + Treatment + Fusion.Partners + Richness + Day:Treatment + Day:Fusion.Partners + Day:Richness + Day:Treatment:Fusion.Partners + Day:Treatment:Richness + Day:Fusion.Partners:Richness:Treatment + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family=binomial, subset=c(Community=="Multiple"))
#fixed_form<-glm(Survivorship ~ Day * Treatment * Fusion.Partners *Richness, data=stress, family=binomial, subset=c(Community=="Multiple"))
#library("brglm2"); glm(fixed_form, data = stress, family = binomial, method="detect_separation")
summary(survmod3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Survivorship ~ Day + Treatment + Fusion.Partners + Richness +
## Day:Treatment + Day:Fusion.Partners + Day:Richness + Day:Treatment:Fusion.Partners +
## Day:Treatment:Richness + Day:Fusion.Partners:Richness:Treatment +
## (1 | Parent.Site/Colony) + (1 | Treatment:Tank) + (1 | Slide/Sample)
## Data: stress
## Subset: c(Community == "Multiple")
##
## AIC BIC logLik deviance df.resid
## 222.5 301.5 -94.2 188.5 756
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.6798 0.0000 0.0000 0.0000 0.5145
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample:Slide (Intercept) 940.14567 30.6618
## Slide (Intercept) 0.04049 0.2012
## Colony:Parent.Site (Intercept) 0.01752 0.1324
## Treatment:Tank (Intercept) 8.33415 2.8869
## Parent.Site (Intercept) 0.03312 0.1820
## Number of obs: 773, groups:
## Sample:Slide, 111; Slide, 45; Colony:Parent.Site, 20; Treatment:Tank, 8; Parent.Site, 6
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) -2.29443 39.97253 -0.057
## Day 0.55107 1.24309 0.443
## TreatmentHigh 111.46301 35.69621 3.123
## Fusion.Partners 100.44513 51.46355 1.952
## Richness 22.39796 34.13494 0.656
## Day:TreatmentHigh -4.98049 1.87839 -2.651
## Day:Fusion.Partners -2.73031 1.88657 -1.447
## Day:Richness -0.90135 1.21483 -0.742
## Day:TreatmentHigh:Fusion.Partners -1.42109 0.99729 -1.425
## Day:TreatmentHigh:Richness -0.02063 0.80677 -0.026
## Day:TreatmentAmbient:Fusion.Partners:Richness -0.09927 0.27596 -0.360
## Day:TreatmentHigh:Fusion.Partners:Richness 0.07110 0.22617 0.314
## Pr(>|z|)
## (Intercept) 0.95423
## Day 0.65754
## TreatmentHigh 0.00179 **
## Fusion.Partners 0.05097 .
## Richness 0.51172
## Day:TreatmentHigh 0.00801 **
## Day:Fusion.Partners 0.14783
## Day:Richness 0.45811
## Day:TreatmentHigh:Fusion.Partners 0.15417
## Day:TreatmentHigh:Richness 0.97959
## Day:TreatmentAmbient:Fusion.Partners:Richness 0.71905
## Day:TreatmentHigh:Fusion.Partners:Richness 0.75326
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Day TrtmnH Fsn.Pr Rchnss Dy:TrH Dy:F.P Dy:Rch Dy:TH:F.P
## Day -0.633
## TreatmntHgh 0.030 -0.317
## Fusn.Prtnrs 0.516 -0.509 0.005
## Richness -0.917 0.826 -0.219 -0.590
## Dy:TrtmntHg -0.392 0.059 -0.623 -0.043 0.331
## Dy:Fsn.Prtn -0.611 0.393 0.087 -0.936 0.590 0.133
## Day:Richnss 0.627 -0.990 0.337 0.519 -0.838 -0.074 -0.404
## Dy:TrtH:F.P 0.061 0.336 -0.155 -0.324 0.137 -0.338 0.071 -0.337
## Dy:TrtmnH:R 0.478 0.213 -0.071 0.098 -0.287 -0.680 -0.275 -0.210 0.453
## Dy:TA:F.P:R 0.424 0.315 -0.245 0.088 -0.181 -0.371 -0.391 -0.313 0.552
## Dy:TH:F.P:R 0.162 -0.148 -0.066 0.147 -0.178 0.307 -0.142 0.151 -0.652
## D:TH:R D:TA:F
## Day
## TreatmntHgh
## Fusn.Prtnrs
## Richness
## Dy:TrtmntHg
## Dy:Fsn.Prtn
## Day:Richnss
## Dy:TrtH:F.P
## Dy:TrtmnH:R
## Dy:TA:F.P:R 0.752
## Dy:TH:F.P:R -0.255 0.031
## convergence code: 0
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 6 negative eigenvalues
Anova(survmod3, type="II")
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Survivorship
## Chisq Df Pr(>Chisq)
## Day 8.3067 1 0.00395 **
## Treatment 5.5318 1 0.01867 *
## Fusion.Partners 0.6676 1 0.41390
## Richness 0.2974 1 0.58553
## Day:Treatment 16.8143 1 4.122e-05 ***
## Day:Fusion.Partners 0.8008 1 0.37084
## Day:Richness 1.0944 1 0.29549
## Day:Treatment:Fusion.Partners 4.1122 1 0.04257 *
## Day:Treatment:Richness 0.3144 1 0.57499
## Day:Treatment:Fusion.Partners:Richness 0.2355 2 0.88894
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dispersion_glmer(survmod3) #no overdispersion
## [1] 0.2731621
plot(residuals(survmod3)) + abline(h=0, lty=2) #Dispersed randomly
## integer(0)
meanstresssurv2 <- ddply(final_stress, c("Fusion.Partners", "Treatment"), summarise,
N = length(Survivorship[!is.na(Survivorship)]),
mean = mean(Survivorship, na.rm=TRUE),
sd = sd(Survivorship, na.rm=TRUE),
se = sd / sqrt(N),
max = max(Survivorship, na.rm=TRUE),
lower = mean-sd,
upper = mean+sd
)
meanstresssurv2
## Fusion.Partners Treatment N mean sd se max lower
## 1 0 Ambient 30 0.83333333 0.3790490 0.06920457 1 0.4542843
## 2 0 High 50 0.04000000 0.1979487 0.02799417 1 -0.1579487
## 3 1 Ambient 16 0.81250000 0.4031129 0.10077822 1 0.4093871
## 4 1 High 15 0.06666667 0.2581989 0.06666667 1 -0.1915322
## 5 2 Ambient 9 0.77777778 0.4409586 0.14698618 1 0.3368192
## 6 2 High 16 0.00000000 0.0000000 0.00000000 0 0.0000000
## 7 3 Ambient 12 0.83333333 0.3892495 0.11236664 1 0.4440839
## 8 3 High 8 0.00000000 0.0000000 0.00000000 0 0.0000000
## upper
## 1 1.2123824
## 2 0.2379487
## 3 1.2156129
## 4 0.3248656
## 5 1.2187363
## 6 0.0000000
## 7 1.2225828
## 8 0.0000000
Plot effect of number of juveniles over time.
stress$group<-paste(stress$Treatment, "-", stress$Community)
survmod3b2<-glmer(Survivorship ~ group * Day + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=stress, family=binomial, link=logit)
eff.surv3b2 <- predictorEffect(c("Day"), survmod3b2)
surv.effects3b2<-plot(eff.surv3b2,
lines=list(multiline=TRUE,
col=c("blue", "blue", "red", "red"),
lty=c(1, 2, 1, 2)),
confint=list(style="bands", alpha=0.1),
lwd=4,
axes=list(y=list(lim=c(0, 1.1))),
type="response",
ylab=expression(bold("Prob(Survivorship)")),
legend.position="top",
xlab=expression(bold("Day")),
main="",
lattice=list(key.args=list(space="top",
border=FALSE,
title=expression(bold("Temperature - Number of Juveniles")),
cex=1,
cex.title=1)));surv.effects3b2
Examine curve estimates to identify 50% survivorship.
summary(effect(c("group", "Day"), xlevels=50, survmod3b2))
##
## group*Day effect
## Day
## group 0 0.653 1.31 1.96 2.61
## Ambient - Individual 0.9989505 0.9988891 0.9988237 0.9987553 0.9986828
## Ambient - Multiple 0.9999639 0.9999570 0.9999487 0.9999389 0.9999273
## High - Individual 0.9995124 0.9993782 0.9992059 0.9989885 0.9987118
## High - Multiple 0.9999991 0.9999986 0.9999980 0.9999972 0.9999960
## Day
## group 3.27 3.92 4.57 5.22 5.88
## Ambient - Individual 0.9986050 0.9985238 0.9984379 0.9983470 0.9982494
## Ambient - Multiple 0.9999132 0.9998966 0.9998769 0.9998534 0.9998250
## High - Individual 0.9983533 0.9979030 0.9973299 0.9966008 0.9956574
## High - Multiple 0.9999942 0.9999916 0.9999880 0.9999828 0.9999752
## Day
## group 6.53 7.18 7.84 8.49 9.14
## Ambient - Individual 0.9981475 0.9980398 0.9979240 0.9978033 0.9976756
## Ambient - Multiple 0.9997916 0.9997518 0.9997037 0.9996472 0.9995799
## High - Individual 0.9944741 0.9929706 0.9910287 0.9885985 0.9855198
## High - Multiple 0.9999644 0.9999489 0.9999263 0.9998942 0.9998482
## Day
## group 9.8 10.4 11.1 11.8 12.4
## Ambient - Individual 0.9975384 0.9974066 0.9972439 0.9970711 0.9969144
## Ambient - Multiple 0.9994984 0.9994107 0.9992889 0.9991419 0.9989919
## High - Individual 0.9815577 0.9770442 0.9704068 0.9619250 0.9528374
## High - Multiple 0.9997810 0.9996944 0.9995493 0.9993353 0.9990726
## Day
## group 13.1 13.7 14.4 15 15.7
## Ambient - Individual 0.9967210 0.9965456 0.9963292 0.9961329 0.9958907
## Ambient - Multiple 0.9987836 0.9985711 0.9982760 0.9979751 0.9975572
## High - Individual 0.9396333 0.9256365 0.9055717 0.8846468 0.8552516
## High - Multiple 0.9986325 0.9980927 0.9971889 0.9960815 0.9942302
## Day
## group 16.3 17 17.6 18.3 18.9
## Ambient - Individual 0.9956711 0.9954001 0.9951544 0.9948512 0.9945764
## Ambient - Multiple 0.9971311 0.9965396 0.9959368 0.9951002 0.9942481
## High - Individual 0.8253269 0.7844978 0.7443186 0.6916290 0.6420357
## High - Multiple 0.9919665 0.9881940 0.9836013 0.9759955 0.9668167
## Day
## group 19.6 20.2 20.9 21.6 22.2
## Ambient - Individual 0.9942373 0.9939299 0.9935506 0.9931477 0.9927826
## Ambient - Multiple 0.9930663 0.9918634 0.9901965 0.9881921 0.9861554
## High - Individual 0.5801574 0.5249503 0.4598598 0.3961105 0.3440645
## High - Multiple 0.9518077 0.9340054 0.9056046 0.8667256 0.8233276
## Day
## group 22.9 23.5 24.2 24.8 25.5
## Ambient - Individual 0.9923322 0.9919240 0.9914204 0.9909641 0.9904013
## Ambient - Multiple 0.9833386 0.9804813 0.9765376 0.9725465 0.9670535
## High - Individual 0.2878140 0.2442424 0.1993517 0.1660498 0.1330014
## High - Multiple 0.7595582 0.6935999 0.6054450 0.5237196 0.4270620
## Day
## group 26.1 26.8 27.4 28.1 28.7
## Ambient - Individual 0.9898913 0.98926241 0.98869264 0.98799004 0.9873536
## Ambient - Multiple 0.9615127 0.95391658 0.94628912 0.93588877 0.9255105
## High - Individual 0.1092707 0.08635276 0.07027067 0.05502712 0.0444949
## High - Multiple 0.3481675 0.26582641 0.20600815 0.14957280 0.1119266
## Day
## group 29.4 30 30.7 31.3 32
## Ambient - Individual 0.98656895 0.98585830 0.98498226 0.98418900 0.98321132
## Ambient - Multiple 0.91146307 0.89756352 0.87893523 0.86070948 0.83659930
## High - Individual 0.03463450 0.02789024 0.02162627 0.01736948 0.01343575
## High - Multiple 0.07870967 0.05768922 0.03984634 0.02887951 0.01976041
##
## Lower 95 Percent Confidence Limits
## Day
## group 0 0.653 1.31 1.96 2.61
## Ambient - Individual 0.9688623 0.9690922 0.9693058 0.9694989 0.9696732
## Ambient - Multiple 0.9981187 0.9979334 0.9977283 0.9975049 0.9972591
## High - Individual 0.9920845 0.9906655 0.9889803 0.9870133 0.9846947
## High - Multiple 0.9999701 0.9999606 0.9999479 0.9999313 0.9999094
## Day
## group 3.27 3.92 4.57 5.22 5.88
## Ambient - Individual 0.9698303 0.9699644 0.9700772 0.9701676 0.9702353
## Ambient - Multiple 0.9969842 0.9966859 0.9963574 0.9959955 0.9955902
## High - Individual 0.9819168 0.9786896 0.9748889 0.9704148 0.9650636
## High - Multiple 0.9998801 0.9998418 0.9997913 0.9997246 0.9996350
## Day
## group 6.53 7.18 7.84 8.49 9.14
## Ambient - Individual 0.9702769 0.9702921 0.9702789 0.9702358 0.9701607
## Ambient - Multiple 0.9951498 0.9946641 0.9941193 0.9935266 0.9928718
## High - Individual 0.9588613 0.9515783 0.9428952 0.9328710 0.9211579
## High - Multiple 0.9995181 0.9993637 0.9991558 0.9988845 0.9985255
## Day
## group 9.8 10.4 11.1 11.8 12.4
## Ambient - Individual 0.9700495 0.9699156 0.9697168 0.9694682 0.9692122
## Ambient - Multiple 0.9921364 0.9913992 0.9904473 0.9893846 0.9883749
## High - Individual 0.9072779 0.8926825 0.8729571 0.8499648 0.8273689
## High - Multiple 0.9980419 0.9974650 0.9965725 0.9953636 0.9939910
## Day
## group 13.1 13.7 14.4 15 15.7
## Ambient - Individual 0.9688590 0.9685055 0.9680282 0.9675585 0.9669328
## Ambient - Multiple 0.9870676 0.9858233 0.9842088 0.9826685 0.9806649
## High - Individual 0.7973226 0.7681915 0.7300762 0.6938098 0.6473968
## High - Multiple 0.9918652 0.9894505 0.9857101 0.9814630 0.9748907
## Day
## group 16.3 17 17.6 18.3 18.9
## Ambient - Individual 0.9663236 0.9655190 0.9647409 0.9637193 0.9627359
## Ambient - Multiple 0.9787481 0.9762472 0.9738466 0.9707029 0.9676733
## High - Individual 0.6043419 0.5508357 0.5028190 0.4453401 0.3958471
## High - Multiple 0.9674434 0.9559596 0.9430171 0.9232269 0.9011890
## Day
## group 19.6 20.2 20.9 21.6 22.2
## Ambient - Individual 0.9614496 0.9602154 0.9586051 0.9567911 0.9550561
## Ambient - Multiple 0.9636884 0.9598297 0.9547280 0.9488733 0.9431543
## High - Individual 0.3392213 0.2927669 0.2422482 0.1967494 0.1622952
## High - Multiple 0.8680761 0.8320749 0.7797593 0.7153722 0.6506739
## Day
## group 22.9 23.5 24.2 24.8 25.5
## Ambient - Individual 0.9527987 0.9506428 0.94784145 0.94516956 0.9417021
## Ambient - Multiple 0.9355228 0.9280213 0.91794801 0.90798712 0.8945388
## High - Individual 0.1276564 0.1026581 0.07860727 0.06193355 0.0464441
## High - Multiple 0.5657819 0.4877617 0.39554313 0.31998675 0.2407995
## Day
## group 26.1 26.8 27.4 28.1 28.7
## Ambient - Individual 0.93839918 0.93411856 0.93004678 0.92477792 0.91977455
## Ambient - Multiple 0.88118186 0.86309627 0.84511815 0.82082232 0.79679637
## High - Individual 0.03603222 0.02660821 0.02041373 0.01490956 0.01134799
## High - Multiple 0.18322882 0.12923146 0.09371102 0.06306962 0.04427288
## Day
## group 29.4 30 30.7 31.3
## Ambient - Individual 0.913312333 0.907188427 0.899297391 0.891838810
## Ambient - Multiple 0.764632262 0.733264999 0.692078154 0.652889335
## High - Individual 0.008223676 0.006223925 0.004485222 0.003380753
## High - Multiple 0.028915333 0.019893620 0.012757874 0.008671295
## Day
## group 32
## Ambient - Individual 0.882256056
## Ambient - Multiple 0.602984743
## High - Individual 0.002426434
## High - Multiple 0.005499545
##
## Upper 95 Percent Confidence Limits
## Day
## group 0 0.653 1.31 1.96 2.61
## Ambient - Individual 0.9999657 0.9999612 0.9999562 0.9999506 0.9999444
## Ambient - Multiple 0.9999993 0.9999991 0.9999988 0.9999985 0.9999981
## High - Individual 0.9999702 0.9999589 0.9999433 0.9999221 0.9998930
## High - Multiple 1.0000000 1.0000000 0.9999999 0.9999999 0.9999998
## Day
## group 3.27 3.92 4.57 5.22 5.88
## Ambient - Individual 0.9999373 0.9999294 0.9999207 0.9999109 0.9998998
## Ambient - Multiple 0.9999975 0.9999968 0.9999959 0.9999947 0.9999931
## High - Individual 0.9998523 0.9997972 0.9997218 0.9996186 0.9994748
## High - Multiple 0.9999997 0.9999996 0.9999993 0.9999989 0.9999983
## Day
## group 6.53 7.18 7.84 8.49 9.14
## Ambient - Individual 0.9998876 0.9998740 0.9998587 0.9998420 0.9998236
## Ambient - Multiple 0.9999911 0.9999885 0.9999851 0.9999809 0.9999754
## High - Individual 0.9992809 0.9990161 0.9986487 0.9981550 0.9974840
## High - Multiple 0.9999974 0.9999959 0.9999936 0.9999900 0.9999844
## Day
## group 9.8 10.4 11.1 11.8 12.4
## Ambient - Individual 0.9998028 0.9997821 0.9997555 0.9997261 0.9996985
## Ambient - Multiple 0.9999682 0.9999599 0.9999475 0.9999313 0.9999134
## High - Individual 0.9965576 0.9954292 0.9936503 0.9912023 0.9883944
## High - Multiple 0.9999755 0.9999632 0.9999409 0.9999050 0.9998575
## Day
## group 13.1 13.7 14.4 15 15.7
## Ambient - Individual 0.9996634 0.9996306 0.9995892 0.9995507 0.9995024
## Ambient - Multiple 0.9998868 0.9998576 0.9998142 0.9997666 0.9996959
## High - Individual 0.9840224 0.9790595 0.9714308 0.9629020 0.9500353
## High - Multiple 0.9997714 0.9996576 0.9994521 0.9991813 0.9986941
## Day
## group 16.3 17 17.6 18.3 18.9
## Ambient - Individual 0.9994579 0.9994024 0.9993517 0.9992890 0.9992323
## Ambient - Multiple 0.9996189 0.9995047 0.9993806 0.9991973 0.9989992
## High - Individual 0.9359643 0.9152970 0.8933854 0.8623573 0.8307878
## High - Multiple 0.9980548 0.9969114 0.9954211 0.9927783 0.9893703
## Day
## group 19.6 20.2 20.9 21.6 22.2
## Ambient - Individual 0.9991628 0.9991006 0.9990252 0.9989470 0.9988782
## Ambient - Multiple 0.9987079 0.9983946 0.9979371 0.9973571 0.9967406
## High - Individual 0.7881158 0.7468255 0.6939347 0.6372250 0.5868072
## High - Multiple 0.9834108 0.9758592 0.9629580 0.9439063 0.9210072
## Day
## group 22.9 23.5 24.2 24.8 25.5
## Ambient - Individual 0.9987962 0.9987249 0.9986410 0.9985688 0.9984850
## Ambient - Multiple 0.9958518 0.9949165 0.9935835 0.9921980 0.9902508
## High - Individual 0.5274205 0.4772442 0.4208513 0.3751904 0.3257642
## High - Multiple 0.8845106 0.8432974 0.7825316 0.7198529 0.6365925
## Day
## group 26.1 26.8 27.4 28.1 28.7
## Ambient - Individual 0.9984139 0.9983324 0.9982640 0.9981867 0.9981227
## Ambient - Multiple 0.9882570 0.9855000 0.9827247 0.9789551 0.9752286
## High - Individual 0.2870447 0.2463006 0.2151503 0.1830334 0.1589008
## High - Multiple 0.5598160 0.4690315 0.3943239 0.3148500 0.2553418
## Day
## group 29.4 30 30.7 31.3 32
## Ambient - Individual 0.9980511 0.9979928 0.99792836 0.99787649 0.99782004
## Ambient - Multiple 0.9702584 0.9654314 0.95910208 0.95305189 0.94523397
## High - Individual 0.1343732 0.1161636 0.09783686 0.08434176 0.07084932
## High - Multiple 0.1968694 0.1558722 0.11759973 0.09182032 0.06845565
Create plot for effect of temperature on survival in multiple juvenile groups across number of fusion partners.
stress$group<-paste(stress$Treatment, "-", stress$Fusion.Partners)
survmod3a<-glmer(Survivorship ~ Day * group + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family=binomial, link=logit)
summary(effect(c("Day", "group"), xlevels=50, survmod3a))
##
## Day*group effect
## group
## Day Ambient - 0 Ambient - 1 Ambient - 2 Ambient - 3 High - 0
## 0 0.9999998 1.0000000 1.0000000 1.0000000 0.999999361
## 0.653 0.9999998 1.0000000 1.0000000 1.0000000 0.999999040
## 1.31 0.9999998 1.0000000 1.0000000 1.0000000 0.999998555
## 1.96 0.9999997 1.0000000 1.0000000 1.0000000 0.999997835
## 2.61 0.9999996 1.0000000 1.0000000 1.0000000 0.999996754
## 3.27 0.9999995 1.0000000 1.0000000 1.0000000 0.999995104
## 3.92 0.9999994 1.0000000 1.0000000 1.0000000 0.999992661
## 4.57 0.9999992 1.0000000 0.9999999 1.0000000 0.999988998
## 5.22 0.9999990 1.0000000 0.9999999 1.0000000 0.999983508
## 5.88 0.9999988 0.9999999 0.9999999 1.0000000 0.999975124
## 6.53 0.9999985 0.9999999 0.9999999 1.0000000 0.999962712
## 7.18 0.9999981 0.9999999 0.9999998 1.0000000 0.999944106
## 7.84 0.9999976 0.9999999 0.9999997 1.0000000 0.999915693
## 8.49 0.9999969 0.9999998 0.9999997 1.0000000 0.999873628
## 9.14 0.9999961 0.9999997 0.9999995 1.0000000 0.999810579
## 9.8 0.9999951 0.9999996 0.9999994 1.0000000 0.999714310
## 10.4 0.9999940 0.9999995 0.9999992 1.0000000 0.999584931
## 11.1 0.9999923 0.9999993 0.9999988 1.0000000 0.999358274
## 11.8 0.9999901 0.9999991 0.9999984 1.0000000 0.999007970
## 12.4 0.9999877 0.9999987 0.9999978 1.0000000 0.998559176
## 13.1 0.9999843 0.9999982 0.9999970 1.0000000 0.997773633
## 13.7 0.9999806 0.9999976 0.9999960 1.0000000 0.996768234
## 14.4 0.9999751 0.9999967 0.9999944 1.0000000 0.995011139
## 15 0.9999692 0.9999956 0.9999925 1.0000000 0.992767277
## 15.7 0.9999605 0.9999939 0.9999896 0.9999999 0.988859195
## 16.3 0.9999511 0.9999919 0.9999862 0.9999999 0.983893136
## 17 0.9999373 0.9999886 0.9999807 0.9999998 0.975309313
## 17.6 0.9999225 0.9999849 0.9999744 0.9999997 0.964520051
## 18.3 0.9999006 0.9999789 0.9999643 0.9999995 0.946176534
## 18.9 0.9998771 0.9999719 0.9999526 0.9999992 0.923653441
## 19.6 0.9998424 0.9999607 0.9999339 0.9999986 0.886664147
## 20.2 0.9998051 0.9999477 0.9999122 0.9999978 0.843360175
## 20.9 0.9997501 0.9999269 0.9998777 0.9999963 0.776867110
## 21.6 0.9996797 0.9998979 0.9998297 0.9999937 0.692441518
## 22.2 0.9996038 0.9998640 0.9997737 0.9999901 0.607756831
## 22.9 0.9994923 0.9998101 0.9996848 0.9999831 0.500487404
## 23.5 0.9993719 0.9997472 0.9995812 0.9999733 0.408126431
## 24.2 0.9991951 0.9996469 0.9994167 0.9999546 0.308389509
## 24.8 0.9990044 0.9995298 0.9992252 0.9999284 0.234814050
## 25.5 0.9987242 0.9993434 0.9989209 0.9998782 0.165581961
## 26.1 0.9984222 0.9991258 0.9985668 0.9998078 0.120158099
## 26.8 0.9979785 0.9987794 0.9980046 0.9996730 0.081146008
## 27.4 0.9975003 0.9983752 0.9973505 0.9994843 0.057294852
## 28.1 0.9967982 0.9977321 0.9963130 0.9991227 0.037815540
## 28.7 0.9960420 0.9969822 0.9951072 0.9986169 0.026335429
## 29.4 0.9949323 0.9957899 0.9931971 0.9976483 0.017189882
## 30 0.9937382 0.9944014 0.9909816 0.9962955 0.011893945
## 30.7 0.9919878 0.9921975 0.9874811 0.9937118 0.007723732
## 31.3 0.9901067 0.9896366 0.9834351 0.9901170 0.005328358
## 32 0.9873542 0.9855842 0.9770732 0.9832961 0.003452106
## group
## Day High - 1 High - 2 High - 3
## 0 1.000000e+00 1.000000e+00 1.000000e+00
## 0.653 1.000000e+00 1.000000e+00 1.000000e+00
## 1.31 1.000000e+00 1.000000e+00 1.000000e+00
## 1.96 1.000000e+00 1.000000e+00 1.000000e+00
## 2.61 1.000000e+00 1.000000e+00 1.000000e+00
## 3.27 1.000000e+00 1.000000e+00 1.000000e+00
## 3.92 1.000000e+00 1.000000e+00 1.000000e+00
## 4.57 1.000000e+00 1.000000e+00 1.000000e+00
## 5.22 1.000000e+00 1.000000e+00 1.000000e+00
## 5.88 1.000000e+00 1.000000e+00 1.000000e+00
## 6.53 1.000000e+00 1.000000e+00 1.000000e+00
## 7.18 1.000000e+00 1.000000e+00 1.000000e+00
## 7.84 1.000000e+00 1.000000e+00 1.000000e+00
## 8.49 1.000000e+00 1.000000e+00 1.000000e+00
## 9.14 1.000000e+00 1.000000e+00 1.000000e+00
## 9.8 1.000000e+00 1.000000e+00 1.000000e+00
## 10.4 1.000000e+00 1.000000e+00 1.000000e+00
## 11.1 1.000000e+00 1.000000e+00 1.000000e+00
## 11.8 1.000000e+00 1.000000e+00 1.000000e+00
## 12.4 1.000000e+00 1.000000e+00 1.000000e+00
## 13.1 1.000000e+00 1.000000e+00 1.000000e+00
## 13.7 1.000000e+00 1.000000e+00 1.000000e+00
## 14.4 1.000000e+00 1.000000e+00 1.000000e+00
## 15 1.000000e+00 1.000000e+00 1.000000e+00
## 15.7 9.999999e-01 1.000000e+00 1.000000e+00
## 16.3 9.999997e-01 1.000000e+00 1.000000e+00
## 17 9.999988e-01 1.000000e+00 1.000000e+00
## 17.6 9.999964e-01 1.000000e+00 9.999999e-01
## 18.3 9.999867e-01 9.999999e-01 9.999997e-01
## 18.9 9.999594e-01 9.999997e-01 9.999987e-01
## 19.6 9.998514e-01 9.999982e-01 9.999938e-01
## 20.2 9.995479e-01 9.999914e-01 9.999758e-01
## 20.9 9.983464e-01 9.999486e-01 9.998823e-01
## 21.6 9.939706e-01 9.996914e-01 9.994275e-01
## 22.2 9.818793e-01 9.985668e-01 9.977812e-01
## 22.9 9.366914e-01 9.914545e-01 9.892933e-01
## 23.5 8.294440e-01 9.614765e-01 9.596799e-01
## 24.2 5.704313e-01 8.060481e-01 8.302386e-01
## 24.8 3.038502e-01 4.720196e-01 5.574822e-01
## 25.5 1.064895e-01 1.295768e-01 2.056284e-01
## 26.1 3.769687e-02 3.103005e-02 6.251173e-02
## 26.8 1.058335e-02 5.304165e-03 1.351589e-02
## 27.4 3.503529e-03 1.145788e-03 3.516907e-03
## 28.1 9.590995e-04 1.909737e-04 7.246614e-04
## 28.7 3.154492e-04 4.108791e-05 1.867692e-04
## 29.4 8.615482e-05 6.842007e-06 3.838209e-05
## 30 2.831984e-05 1.471842e-06 9.887291e-06
## 30.7 7.733039e-06 2.450847e-07 2.031607e-06
## 31.3 2.541784e-06 5.272197e-08 5.233313e-07
## 32 6.940485e-07 8.779024e-09 1.075316e-07
##
## Lower 95 Percent Confidence Limits
## group
## Day Ambient - 0 Ambient - 1 Ambient - 2 Ambient - 3 High - 0
## 0 0.7502997 0.8837796 0.2726660 0.001463071 9.781682e-01
## 0.653 0.7530553 0.8850246 0.2850875 0.001755680 9.756751e-01
## 1.31 0.7557811 0.8862363 0.2978390 0.002108716 9.728836e-01
## 1.96 0.7584307 0.8873939 0.3106826 0.002527219 9.698123e-01
## 2.61 0.7610327 0.8885096 0.3237284 0.003028014 9.663998e-01
## 3.27 0.7636249 0.8895984 0.3371546 0.003637109 9.625488e-01
## 3.92 0.7661278 0.8906264 0.3505255 0.004355243 9.583357e-01
## 4.57 0.7685800 0.8916089 0.3640133 0.005213359 9.536620e-01
## 5.22 0.7709802 0.8925446 0.3775862 0.006238124 9.484806e-01
## 5.88 0.7733628 0.8934451 0.3914209 0.007481614 9.426474e-01
## 6.53 0.7756543 0.8942812 0.4050634 0.008943861 9.362830e-01
## 7.18 0.7778898 0.8950650 0.4186876 0.010685947 9.292435e-01
## 7.84 0.7801005 0.8958049 0.4324652 0.012794095 9.213387e-01
## 8.49 0.7822178 0.8964758 0.4459413 0.015265370 9.127380e-01
## 9.14 0.7842735 0.8970865 0.4592876 0.018199108 9.032534e-01
## 9.8 0.7862959 0.8976415 0.4726674 0.021734659 8.926374e-01
## 10.4 0.7880753 0.8980858 0.4846472 0.025517218 8.820515e-01
## 11.1 0.7900776 0.8985272 0.4983610 0.030729717 8.684797e-01
## 11.8 0.7919971 0.8988798 0.5117462 0.036947118 8.534803e-01
## 12.4 0.7935734 0.8991062 0.5229211 0.043204447 8.393930e-01
## 13.1 0.7953280 0.8992747 0.5355677 0.051753663 8.214044e-01
## 13.7 0.7967562 0.8993307 0.5460375 0.060299403 8.045510e-01
## 14.4 0.7983290 0.8992835 0.5577763 0.071883203 7.830771e-01
## 15 0.7995927 0.8991378 0.5673935 0.083358947 7.629926e-01
## 15.7 0.8009625 0.8988329 0.5780486 0.098754114 7.374256e-01
## 16.3 0.8020412 0.8984439 0.5866577 0.113828785 7.135083e-01
## 17 0.8031808 0.8978243 0.5960390 0.133783734 6.830038e-01
## 17.6 0.8040485 0.8971351 0.6034671 0.153034006 6.543479e-01
## 18.3 0.8049233 0.8961231 0.6113578 0.178090101 6.175250e-01
## 18.9 0.8055462 0.8950548 0.6174021 0.201815457 5.825359e-01
## 19.6 0.8061114 0.8935411 0.6235392 0.232062422 5.368340e-01
## 20.2 0.8064452 0.8919822 0.6279456 0.260063602 4.925175e-01
## 20.9 0.8066414 0.8898105 0.6319903 0.294890679 4.334194e-01
## 21.6 0.8066077 0.8871956 0.6347280 0.331636339 3.651950e-01
## 22.2 0.8063754 0.8845402 0.6359151 0.364245275 2.999984e-01
## 22.9 0.8058392 0.8808703 0.6357809 0.403008271 2.206979e-01
## 23.5 0.8051256 0.8771447 0.6342008 0.436321444 1.570771e-01
## 24.2 0.8039579 0.8719840 0.6304083 0.474595261 9.667846e-02
## 24.8 0.8026323 0.8667220 0.6252517 0.506313359 5.979171e-02
## 25.5 0.8006517 0.8593858 0.6166659 0.541335748 3.226826e-02
## 26.1 0.7985283 0.8518451 0.6067684 0.569083629 1.838364e-02
## 26.8 0.7954743 0.8412302 0.5917764 0.598104957 9.286373e-03
## 27.4 0.7922833 0.8302030 0.5755222 0.619512541 5.092257e-03
## 28.1 0.7877732 0.8145039 0.5519775 0.639614825 2.495599e-03
## 28.7 0.7831143 0.7980134 0.5273587 0.651825176 1.344199e-03
## 29.4 0.7765743 0.7743011 0.4928932 0.658781792 6.490411e-04
## 30 0.7698412 0.7492115 0.4581229 0.656774249 3.463527e-04
## 30.7 0.7603954 0.7130463 0.4114705 0.641900272 1.658714e-04
## 31.3 0.7506579 0.6749697 0.3667811 0.614466495 8.803733e-05
## 32 0.7369611 0.6209891 0.3106283 0.558769702 4.195190e-05
## group
## Day High - 1 High - 2 High - 3
## 0 9.574248e-01 1.000000e+00 9.997589e-01
## 0.653 9.536663e-01 1.000000e+00 9.997065e-01
## 1.31 9.495633e-01 1.000000e+00 9.996422e-01
## 1.96 9.451623e-01 1.000000e+00 9.995647e-01
## 2.61 9.403956e-01 1.000000e+00 9.994704e-01
## 3.27 9.351535e-01 1.000000e+00 9.993537e-01
## 3.92 9.295660e-01 1.000000e+00 9.992134e-01
## 4.57 9.235270e-01 1.000000e+00 9.990426e-01
## 5.22 9.170050e-01 9.999999e-01 9.988345e-01
## 5.88 9.098551e-01 9.999999e-01 9.985766e-01
## 6.53 9.022596e-01 9.999998e-01 9.982666e-01
## 7.18 8.940795e-01 9.999997e-01 9.978886e-01
## 7.84 8.851382e-01 9.999995e-01 9.974197e-01
## 8.49 8.756696e-01 9.999991e-01 9.968555e-01
## 9.14 8.655056e-01 9.999984e-01 9.961670e-01
## 9.8 8.544346e-01 9.999973e-01 9.953120e-01
## 10.4 8.436800e-01 9.999955e-01 9.943687e-01
## 11.1 8.302604e-01 9.999920e-01 9.930232e-01
## 11.8 8.158546e-01 9.999856e-01 9.913524e-01
## 12.4 8.026845e-01 9.999763e-01 9.896011e-01
## 13.1 7.863152e-01 9.999576e-01 9.870983e-01
## 13.7 7.713860e-01 9.999301e-01 9.844711e-01
## 14.4 7.528711e-01 9.998746e-01 9.807097e-01
## 15 7.360161e-01 9.997929e-01 9.767532e-01
## 15.7 7.151396e-01 9.996277e-01 9.710741e-01
## 16.3 6.961435e-01 9.993838e-01 9.650815e-01
## 17 6.725964e-01 9.988889e-01 9.564440e-01
## 17.6 6.511163e-01 9.981553e-01 9.472801e-01
## 18.3 6.243514e-01 9.966592e-01 9.339720e-01
## 18.9 5.997180e-01 9.944276e-01 9.197058e-01
## 19.6 5.685690e-01 9.898405e-01 8.986694e-01
## 20.2 5.392395e-01 9.829333e-01 8.756143e-01
## 20.9 5.008013e-01 9.685580e-01 8.404337e-01
## 21.6 4.556198e-01 9.415950e-01 7.918036e-01
## 22.2 4.083691e-01 8.997846e-01 7.327170e-01
## 22.9 3.364591e-01 8.088340e-01 6.275359e-01
## 23.5 2.519680e-01 6.617411e-01 4.827427e-01
## 24.2 1.276788e-01 3.601855e-01 2.394239e-01
## 24.8 4.243118e-02 1.029092e-01 6.943812e-02
## 25.5 6.571658e-03 1.004368e-02 7.999050e-03
## 26.1 1.013236e-03 9.820130e-04 9.023227e-04
## 26.8 9.959501e-05 5.686141e-05 5.962968e-05
## 27.4 1.286033e-05 4.691155e-06 5.391782e-06
## 28.1 1.138322e-06 2.473758e-07 3.114505e-07
## 28.7 1.398032e-07 1.954346e-08 2.637118e-08
## 29.4 1.193863e-08 9.996445e-10 1.452553e-09
## 30 1.437237e-09 7.766146e-11 1.197516e-10
## 30.7 1.207823e-10 3.919710e-12 6.455679e-12
## 31.3 1.439865e-11 3.020606e-13 5.252472e-13
## 32 1.199847e-12 2.220446e-16 2.220446e-16
##
## Upper 95 Percent Confidence Limits
## group
## Day Ambient - 0 Ambient - 1 Ambient - 2 Ambient - 3 High - 0 High - 1
## 0 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 0.653 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 1.31 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 1.96 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 2.61 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 3.27 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 3.92 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 4.57 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 5.22 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 5.88 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 6.53 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 7.18 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## 7.84 1.0000000 1.0000000 1.0000000 1.0000000 0.9999999 1.0000000
## 8.49 1.0000000 1.0000000 1.0000000 1.0000000 0.9999998 1.0000000
## 9.14 1.0000000 1.0000000 1.0000000 1.0000000 0.9999997 1.0000000
## 9.8 1.0000000 1.0000000 1.0000000 1.0000000 0.9999993 1.0000000
## 10.4 1.0000000 1.0000000 1.0000000 1.0000000 0.9999987 1.0000000
## 11.1 1.0000000 1.0000000 1.0000000 1.0000000 0.9999973 1.0000000
## 11.8 1.0000000 1.0000000 1.0000000 1.0000000 0.9999943 1.0000000
## 12.4 1.0000000 1.0000000 1.0000000 1.0000000 0.9999891 1.0000000
## 13.1 1.0000000 1.0000000 1.0000000 1.0000000 0.9999771 1.0000000
## 13.7 1.0000000 1.0000000 1.0000000 1.0000000 0.9999567 1.0000000
## 14.4 1.0000000 1.0000000 1.0000000 1.0000000 0.9999093 1.0000000
## 15 1.0000000 1.0000000 1.0000000 1.0000000 0.9998292 1.0000000
## 15.7 1.0000000 1.0000000 1.0000000 1.0000000 0.9996437 1.0000000
## 16.3 1.0000000 1.0000000 1.0000000 1.0000000 0.9993330 1.0000000
## 17 1.0000000 1.0000000 1.0000000 1.0000000 0.9986210 1.0000000
## 17.6 1.0000000 1.0000000 1.0000000 1.0000000 0.9974449 1.0000000
## 18.3 1.0000000 1.0000000 1.0000000 1.0000000 0.9948026 1.0000000
## 18.9 0.9999999 1.0000000 1.0000000 1.0000000 0.9905563 1.0000000
## 19.6 0.9999999 1.0000000 1.0000000 1.0000000 0.9814146 1.0000000
## 20.2 0.9999998 1.0000000 1.0000000 1.0000000 0.9676051 0.9999998
## 20.9 0.9999997 1.0000000 1.0000000 1.0000000 0.9406387 0.9999972
## 21.6 0.9999996 0.9999999 0.9999999 1.0000000 0.8980738 0.9999692
## 22.2 0.9999993 0.9999999 0.9999999 1.0000000 0.8485266 0.9997650
## 22.9 0.9999989 0.9999997 0.9999998 1.0000000 0.7799720 0.9976890
## 23.5 0.9999984 0.9999995 0.9999997 1.0000000 0.7184340 0.9859575
## 24.2 0.9999973 0.9999992 0.9999994 1.0000000 0.6500751 0.9233573
## 24.8 0.9999960 0.9999986 0.9999990 1.0000000 0.5969047 0.8112961
## 25.5 0.9999934 0.9999974 0.9999981 1.0000000 0.5414880 0.6822590
## 26.1 0.9999901 0.9999956 0.9999968 1.0000000 0.4989679 0.6020676
## 26.8 0.9999840 0.9999921 0.9999942 0.9999998 0.4541599 0.5346024
## 27.4 0.9999760 0.9999871 0.9999904 0.9999996 0.4191760 0.4901014
## 28.1 0.9999617 0.9999773 0.9999831 0.9999986 0.3817223 0.4474057
## 28.7 0.9999430 0.9999638 0.9999730 0.9999964 0.3521295 0.4159638
## 29.4 0.9999098 0.9999387 0.9999544 0.9999893 0.3202061 0.3834167
## 30 0.9998672 0.9999053 0.9999300 0.9999735 0.2948768 0.3581746
## 30.7 0.9997930 0.9998464 0.9998876 0.9999282 0.2675132 0.3311540
## 31.3 0.9996995 0.9997723 0.9998357 0.9998412 0.2458115 0.3097267
## 32 0.9995406 0.9996496 0.9997520 0.9996347 0.2224093 0.2864642
## group
## Day High - 2 High - 3
## 0 1.000000000 1.0000000
## 0.653 1.000000000 1.0000000
## 1.31 1.000000000 1.0000000
## 1.96 1.000000000 1.0000000
## 2.61 1.000000000 1.0000000
## 3.27 1.000000000 1.0000000
## 3.92 1.000000000 1.0000000
## 4.57 1.000000000 1.0000000
## 5.22 1.000000000 1.0000000
## 5.88 1.000000000 1.0000000
## 6.53 1.000000000 1.0000000
## 7.18 1.000000000 1.0000000
## 7.84 1.000000000 1.0000000
## 8.49 1.000000000 1.0000000
## 9.14 1.000000000 1.0000000
## 9.8 1.000000000 1.0000000
## 10.4 1.000000000 1.0000000
## 11.1 1.000000000 1.0000000
## 11.8 1.000000000 1.0000000
## 12.4 1.000000000 1.0000000
## 13.1 1.000000000 1.0000000
## 13.7 1.000000000 1.0000000
## 14.4 1.000000000 1.0000000
## 15 1.000000000 1.0000000
## 15.7 1.000000000 1.0000000
## 16.3 1.000000000 1.0000000
## 17 1.000000000 1.0000000
## 17.6 1.000000000 1.0000000
## 18.3 1.000000000 1.0000000
## 18.9 1.000000000 1.0000000
## 19.6 1.000000000 1.0000000
## 20.2 0.999999996 1.0000000
## 20.9 0.999999919 0.9999999
## 21.6 0.999998463 0.9999988
## 22.2 0.999981506 0.9999864
## 22.9 0.999685777 0.9998027
## 23.5 0.996869221 0.9983553
## 24.2 0.968434791 0.9870097
## 24.8 0.874487489 0.9550946
## 25.5 0.685961381 0.8925867
## 26.1 0.510590514 0.8311691
## 26.8 0.333354347 0.7589142
## 27.4 0.219050869 0.6979023
## 28.1 0.128530983 0.6280506
## 28.7 0.079520046 0.5695691
## 29.4 0.044735372 0.5035452
## 30 0.027137480 0.4494486
## 30.7 0.015092945 0.3900025
## 31.3 0.009118243 0.3427206
## 32 0.005064860 0.2923095
eff.surv3 <- predictorEffect(c("Day"), survmod3a)
surv.effects3<-plot(eff.surv3,
lines=list(multiline=TRUE,
col=c("blue", "blue", "blue", "blue", "red", "red", "red", "red"),
lty=c(1, 2, 4, 3, 1, 2, 4, 3)),
confint=list(style="bands", alpha=0.1),
lwd=4,
axes=list(y=list(lim=c(0, 1.1))),
type="response",
ylab=expression(bold("Prob(Survivorship)")),
legend.position="top",
xlab=expression(bold("Day")),
main="",
lattice=list(key.args=list(space="top",
border=FALSE,
title=expression(bold("Temperature - Fusion Partners")),
cex=1,
cex.title=1)));surv.effects3
Display mean survivorship within multiple juvenile groups at the end of the experiment.
meanstresssurv2 <- ddply(multiple_final_stress, c("Fusion.Partners", "Treatment"), summarise,
N = length(Survivorship[!is.na(Survivorship)]),
mean = mean(Survivorship, na.rm=TRUE),
sd = sd(Survivorship, na.rm=TRUE),
se = sd / sqrt(N),
max = max(Survivorship, na.rm=TRUE),
lower = mean-sd,
upper = mean+sd
)
meanstresssurv2
## Fusion.Partners Treatment N mean sd se max lower
## 1 0 Ambient 9 0.55555556 0.5270463 0.17568209 1 0.02850928
## 2 0 High 26 0.07692308 0.2717465 0.05329387 1 -0.19482341
## 3 1 Ambient 16 0.81250000 0.4031129 0.10077822 1 0.40938711
## 4 1 High 15 0.06666667 0.2581989 0.06666667 1 -0.19153222
## 5 2 Ambient 9 0.77777778 0.4409586 0.14698618 1 0.33681923
## 6 2 High 16 0.00000000 0.0000000 0.00000000 0 0.00000000
## 7 3 Ambient 12 0.83333333 0.3892495 0.11236664 1 0.44408386
## 8 3 High 8 0.00000000 0.0000000 0.00000000 0 0.00000000
## upper
## 1 1.0826018
## 2 0.3486696
## 3 1.2156129
## 4 0.3248656
## 5 1.2187363
## 6 0.0000000
## 7 1.2225828
## 8 0.0000000
Analyze with a logistic binomial regression model.
Analyze fusion success over the course of the stress test.
#with full
fusemod2<-glmer(Fusion ~ Day * Treatment * Richness + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family=binomial, link=logit, subset=c(Community=="Multiple"))
summary(fusemod2)
Anova(fusemod2, type="II")
plot(allEffects(fusemod2))
dispersion_glmer(fusemod2) #no overdispersion
plot(residuals(fusemod2)) + abline(h=0, lty=2)
meanstressfuse <- ddply(multiple_final_stress, c("Treatment", "Richness", "Timepoint"), summarise,
N = length(Fusion[!is.na(Fusion)]),
mean = mean(Fusion, na.rm=TRUE),
sd = sd(Fusion, na.rm=TRUE),
se = sd / sqrt(N),
max = max(Fusion, na.rm=TRUE),
lower = mean-sd,
upper = mean+sd
)
meanstressfuse
## Treatment Richness Timepoint N mean sd se max
## 1 Ambient 1 End 19 0.7368421 0.4524139 0.10379087 1
## 2 Ambient 2 End 23 0.8260870 0.3875534 0.08081047 1
## 3 Ambient 4 End 4 1.0000000 0.0000000 0.00000000 1
## 4 High 1 End 25 0.6000000 0.5000000 0.10000000 1
## 5 High 2 End 18 0.5555556 0.5113100 0.12051692 1
## 6 High 3 End 14 0.4285714 0.5135526 0.13725270 1
## 7 High 4 End 8 1.0000000 0.0000000 0.00000000 1
## lower upper
## 1 0.28442818 1.189256
## 2 0.43853357 1.213640
## 3 1.00000000 1.000000
## 4 0.10000000 1.100000
## 5 0.04424556 1.066866
## 6 -0.08498116 0.942124
## 7 1.00000000 1.000000
Fusion is decreased at the end of the experiment.
Analyze with a poisson regression model.
Analyze growth during the course of the stress test.
Test for effect of temperature:
#with full dataset
#polypmod2<-glmer(Polyps ~ Treatment + Community + Community:Treatment + (1|Parent.Site/Genotype) + (1|Treatment:Tank) + (1|Slide), data=stress, family = poisson, subset=c(Timepoint=="S2D11"))
#summary(polypmod2)
#Anova(polypmod2, type="II")
#qqPlot(residuals(polypmod2))
polypmod2a<-glmer(Polyps^2 ~ Day * Treatment * Community + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family = poisson)
summary(polypmod2a)
Anova(polypmod2a, type="II")
qqPlot(residuals(polypmod2a))
Analyze within multiple juvenile groups:
polypmod3<-glmer(Polyps^2 ~ Day * Treatment * Fusion.Partners * Richness + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family = poisson, subset=c(Community==c("Multiple")))
summary(polypmod3)
Anova(polypmod3, type="II")
qqPlot(residuals(polypmod3))
Plot difference in temperature between single and multiple groups:
stress$group<-paste(stress$Treatment, "-", stress$Community)
polypmod2b<-glmer(Polyps^2 ~ Day * group + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide/Sample), data=stress, family = poisson)
eff.surv5 <- predictorEffect("Day", xlevels=4, polypmod2b)
polyp.effects2<-plot(eff.surv5,
lines=list(multiline=TRUE,
col=c("blue", "blue", "red", "red"),
lty=c(1,2, 1, 2)),
confint=list(style="bands", width=4),
lwd=4,
#axes=list(y=list(lim=c(0, 1))),
type="response",
ylab=expression(bold("(Total Polyp Expansion)"^2)),
legend.position="right",
xlab=expression(bold("Time (Days)")),
main="",
lattice=list(key.args=list(space="top",
border=FALSE,
title=expression(bold("Treatment - Number of Juveniles")),
cex=1,
cex.title=1)));polyp.effects2
Present means of growth by number of juveniles. Sample size was too low due to mortality at the end of the experiment to calculate growth within multiple juvenile groups, so we pool these together for analysis here.
meanstressgrowth <- ddply(final_stress, c("Community", "Treatment","Fusion.Partners"), summarise,
N = length(Polyps[!is.na(Polyps)]),
mean = mean(Polyps, na.rm=TRUE),
sd = sd(Polyps, na.rm=TRUE),
se = sd / sqrt(N),
max = max(Polyps, na.rm=TRUE),
lower = mean-sd,
upper = mean+sd
)
meanstressgrowth
Combine figures for Stress Period
stress_figs<-plot_grid(surv.effects3b2, surv.effects3, polyp.effects2, labels = c("A", "B", "C"), ncol=3, nrow=1, rel_heights= c(1,1,1), rel_widths = c(1,1,1), label_size = 20, label_y=1, align="h");stress_figs
ggsave(filename="Figures/stress_figs.png", plot=stress_figs, dpi=500, width=20, height=6, units="in")
Load in temperature data from Apex files. Data were recorded every 15 minutes and probes were calibrated to a digital thermometer weekly.
juvtemps<-read.csv("Data/ApexTemps.csv", header=T, na.strings="NA")
#convert Date and Time to date/time format
juvtemps$Date <- as.POSIXct(juvtemps$Date, format="%m/%d/%y %H:%M")
#convert to long format
juvtemps<-juvtemps %>% gather(Tank, Temperature, Tank17:Tank24)
juvtemps$Tank<-as.factor(juvtemps$Tank)
dailymax<-juvtemps[which(juvtemps$TimeDay == "930"),]
dailymax<-juvtemps[which(juvtemps$Phase == "Stress"),]
dailymax$Treatment<-ifelse(dailymax$Tank %in% c("Tank17"), "Ambient",
ifelse(dailymax$Tank %in% c("Tank18"), "Ambient",
ifelse(dailymax$Tank %in% c("Tank19"), "High",
ifelse(dailymax$Tank %in% c("Tank20"), "High",
ifelse(dailymax$Tank %in% c("Tank21"), "Ambient",
ifelse(dailymax$Tank %in% c("Tank22"), "High",
ifelse(dailymax$Tank %in% c("Tank23"), "Ambient",
ifelse(dailymax$Tank %in% c("Tank24"), "High",NA))))))))
maxtemps <- ddply(dailymax, c("Treatment"), summarise,
N = length(Temperature[!is.na(Temperature)]),
mean = mean(Temperature, na.rm=TRUE),
sd = sd(Temperature, na.rm=TRUE),
se = sd / sqrt(N),
max = max(Temperature, na.rm=TRUE),
lower = mean-sd,
upper = mean+sd
)
juvtemps$Period<-ifelse(juvtemps$Phase %in% c("Stress"), "Stress",
ifelse(juvtemps$Phase %in% c("Stress2"), "Stress",
ifelse(juvtemps$Phase %in% c("Acclimation"), "Stress",
ifelse(juvtemps$Phase %in% c("Hurricane"), "Hurricane",
ifelse(juvtemps$Phase %in% c("Ambient"), "Grow-Out",NA)))))
Plot timeseries of temperature data with tank temperature in colors.
#subset temperature measurements
juvtempsHIGH<-subset(juvtemps, Period==c("Stress"), c(Date, TimeDay, Phase, Tank, Temperature, Period))
juvtempsHIGH$Treatment<-ifelse(juvtempsHIGH$Tank %in% c("Tank17"), "Ambient",
ifelse(juvtempsHIGH$Tank %in% c("Tank18"), "Ambient",
ifelse(juvtempsHIGH$Tank %in% c("Tank19"), "High",
ifelse(juvtempsHIGH$Tank %in% c("Tank20"), "High",
ifelse(juvtempsHIGH$Tank %in% c("Tank21"), "Ambient",
ifelse(juvtempsHIGH$Tank %in% c("Tank22"), "High",
ifelse(juvtempsHIGH$Tank %in% c("Tank23"), "Ambient",
ifelse(juvtempsHIGH$Tank %in% c("Tank24"), "High",NA))))))))
juvtempsHURR<-subset(juvtemps, Period==c("Hurricane"), c(Date, TimeDay, Phase, Tank, Temperature, Period))
juvtempsHURR$Treatment<-ifelse(juvtempsHURR$Tank %in% c("Tank17"), "Ambient",
ifelse(juvtempsHURR$Tank %in% c("Tank18"), "Ambient",
ifelse(juvtempsHURR$Tank %in% c("Tank19"), "High",
ifelse(juvtempsHURR$Tank %in% c("Tank20"), "High",
ifelse(juvtempsHURR$Tank %in% c("Tank21"), "Ambient",
ifelse(juvtempsHURR$Tank %in% c("Tank22"), "High",
ifelse(juvtempsHURR$Tank %in% c("Tank23"), "Ambient",
ifelse(juvtempsHURR$Tank %in% c("Tank24"), "High",NA))))))))
juvtempsGROW<-subset(juvtemps, Period==c("Grow-Out"), c(Date, TimeDay, Phase, Tank, Temperature, Period))
juvtempsGROW$Treatment<-ifelse(juvtempsGROW$Tank %in% c("Tank17"), "Ambient",
ifelse(juvtempsGROW$Tank %in% c("Tank18"), "Ambient",
ifelse(juvtempsGROW$Tank %in% c("Tank19"), "Ambient",
ifelse(juvtempsGROW$Tank %in% c("Tank20"), "Ambient",
ifelse(juvtempsGROW$Tank %in% c("Tank21"), "Ambient",
ifelse(juvtempsGROW$Tank %in% c("Tank22"), "Ambient",
ifelse(juvtempsGROW$Tank %in% c("Tank23"), "Ambient",
ifelse(juvtempsGROW$Tank %in% c("Tank24"), "Ambient",NA))))))))
juvtemps<-rbind(juvtempsHIGH, juvtempsGROW)
juvtemps<-rbind(juvtemps, juvtempsHURR)
TimeSeries1j<-ggplot(juvtemps, aes(x = Date, y = Temperature)) +
geom_line(aes(color = Treatment), size = 1) +
scale_color_manual(values=c("blue", "red")) +
ylab("Temperature (°C)")+
xlab("Date")+
ylim(26,32)+
geom_vline(xintercept = as.POSIXct(as.Date(c("2018-07-31 10:00:00"))), linetype="solid",
color = "black", size=1)+
geom_vline(xintercept = as.POSIXct(as.Date(c("2018-08-16 00:00:00"))), linetype="dashed",
color = "black", size=1)+
geom_vline(xintercept = as.POSIXct(as.Date(c("2018-09-20 16:00:00"))), linetype="dotted",
color = "black", size=1)+
theme_classic()+
scale_x_datetime(limits = as.POSIXct(as.Date(c("2018-07-30 09:00:00", "2018-09-20 21:00:00")))) +
theme(legend.position="none")+
theme(text = element_text(size = 18, color="black"))+
theme(axis.text = element_text(size=16, color="black"));TimeSeries1j
Plot by time of day by summarizing the mean temp at each time of day for each tank. Display the mean temperatures of each treatment during the Stress test phase (excluding the pause in treatment due to the hurricane system shutdown).
#subset juv temps for the stress period
juvtempsSTRESS<-juvtemps
juvtempsSTRESS$Period<-ifelse(juvtempsSTRESS$Phase %in% c("Stress"), "Stress",
ifelse(juvtempsSTRESS$Phase %in% c("Stress2"), "Stress",
ifelse(juvtempsSTRESS$Phase %in% c("Acclimation"), "Acclimation",
ifelse(juvtempsSTRESS$Phase %in% c("Hurricane"), "Hurricane",
ifelse(juvtempsSTRESS$Phase %in% c("Ambient"), "Grow-Out",NA)))))
juvtempsSTRESS
juvtempsSTRESS<-juvtempsSTRESS[which(juvtempsSTRESS$Phase == "Stress"),]
QuarterlyJ <- ddply(juvtempsSTRESS, c("TimeDay", "Treatment"), summarise,
N = length(Temperature[!is.na(Temperature)]),
mean = mean(Temperature, na.rm=TRUE),
sd = sd(Temperature, na.rm=TRUE),
se = sd / sqrt(N),
max = max(Temperature, na.rm=TRUE),
lower = mean-sd,
upper = mean+sd
)
MeansJtemp <- ddply(juvtemps, c("Period", "Treatment"), summarise,
N = length(Temperature[!is.na(Temperature)]),
mean = mean(Temperature, na.rm=TRUE),
sd = sd(Temperature, na.rm=TRUE),
se = sd / sqrt(N)
)
MeansJtemp
MeansJtemp2 <- ddply(juvtemps, c("Treatment"), summarise,
N = length(Temperature[!is.na(Temperature)]),
mean = mean(Temperature, na.rm=TRUE),
max = max(Temperature, na.rm=TRUE),
sd = sd(Temperature, na.rm=TRUE),
se = sd / sqrt(N)
)
MeansJtemp2
Generate final figure.
JuvTempsDay<-ggplot(QuarterlyJ, aes(x = TimeDay, y = mean, color=Treatment, fill=Treatment)) +
geom_line(color="black", size = 1) +
#facet_wrap(~Period) +
scale_x_continuous(breaks=c(0, 720, 1440), labels=c("0:00", "12:00", "24:00")) +
scale_color_manual(values = c("blue", "red")) +
scale_fill_manual(values = c("blue", "red")) +
theme_classic() +
geom_ribbon(aes(ymin=QuarterlyJ$lower, ymax=QuarterlyJ$upper), linetype=2, alpha=0.4, color=NA) +
ylab("Temperature (°C)") +
xlab("Time of Day") +
ylim(26,32)+
theme(text = element_text(size = 18, color="black"))+
theme(panel.margin = unit(2, "lines"))+
theme(axis.text = element_text(size=16, color="black")); JuvTempsDay
Combine figures for temperature.
temp_figs<-plot_grid(TimeSeries1j, JuvTempsDay, labels = c("A", "B"), ncol=2, nrow=1, rel_heights= c(1,1), rel_widths = c(1,0.8), label_size = 20, label_y=1, align="h");temp_figs
ggsave(filename="Figures/temp_figs.png", plot=temp_figs, dpi=500, width=14, height=6, units="in")